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ABSTRACT 


This  report  describes  a  detailed  investigation  of  the  effects  of  computational  accuracy  on  the  prediction  of  shock 
wave/boundary  layer  interaction.  In  particular,  the  result  of  inaccuracies  in  the  computation  of  the  shock  and  the 
flow  in  its  vicinity  is  studied.  A  new  computational  procedure  computes  all  shocks  as  discontinuities  while 
including  all  viscous  effects.  This  scheme  is  used  as  a  standard  against  which  the  accuracy  of  widely  used  shock 
capturing  schemes  is  measured.  The  effect  of  the  numerical  error  generated  by  spreading  a  shock  over  a  few  mesh 
intervals  (instead  of  a  few  mean  free  paths)  is  evaluated  with  regard  to  shock/boundary  layer  interaction.  We  consider 
the  spreading  error,  as  well  as  the  error  produced  by  reducing  the  formal  accuracy  of  these  schemes  near  shocks  (in 
order  to  eliminate  wiggles).  In  this  report,  we  present  the  computational  scheme  and  results  for  a  number  of  flow 
configurations. 
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I.  INTRODUCTION 


Over  the  past  two  decades,  computational  fluid  dynamics  (CFD)  has  matured  to  the  point  that  it  contributes 
significantly  to  the  understanding  of  fluid  flows.  Now  our  main  concern  is  the  reliability  of  the  predictions  of  CFD, 
and  therefore  accuracy  of  CFD  algorithms  should  be  the  main  driver  of  our  research.  While  advances  in  computer 
technology  may  enable  CFD  researchers  to  resolve  essentially  smooth  portions  of  complex  flow  fields,  special 
considerations  will  have  to  be  paid  to  "singular  regions"  to  ensure  accuracy.  Adaptive  gridding/grid  clustering  will 
help  resolve  geometric  singularities  (wing  leading  edge,  for  example),  boundary  layers,  vortex  sheets  and  vortex 
centers,  but  to  resolve  shock  waves  may  be  too  much  to  ask  of  any  grid-adapting  scheme.  The  shock  wave  is,  in 
fact,  a  discontinuity  under  the  assumption  of  continuum  flow  (its  thickness  being  on  the  order  of  a  few  mean  free 
paths)  while  boundary  layers,  contact  sheets  and  vortex  centers  spread  as  -yjRtx  with  x  measured  from  their  origin, 

in  laminar  flow  (turbulent  muting  is  faster).  We  offer  the  alternative  of  fitting  the  shocks  in  the  flow  field  while 
resolving  all  other  "singularities”  with  grid  adaptation.  In  this  way  very  high  degrees  of  accuracy  can  be  achieved  for 
the  shock  system  and  the  viscous  regions  as  well  (boundary  layer,  contact  sheet  and  vortices). 

With  this  approach,  we  feel  that  the  accurateAeliable  prediction  of  3-D  unsteady  viscous  flows  is  feasible,  if 
not  yet  with  current  supercomputers,  then  using  the  next  generation  of  multiple  processor  machines.  The  equations 
to  be  solved  are  the  Reynolds-averaged  Navier-Stokes  equations  with  an  appropriate  model  for  the  turbulent  shear 
stresses.  The  computational  procedure  used  in  smooth  portions  of  the  flow  will  assure  consistency  between  the 
numerical  and  physical  domains  of  dependence.  Shock  waves  will  be  fit  and  the  exact  Rankine-Hugoniot  jump 
conditions  satisfied  across  them.  Grid  adaptation  can  be  used  to  resolve  geometric  singularities,  boundary  layers, 
vortex  sheets  and  vortex  centers. 

In  the  phase  of  our  effort  reported  upon  here,  we  were  able  to  demonstrate  the  capability  of  our 
computational  procedure  by  applying  our  scheme  to  a  number  of  unit  problems.  The  ability  to  combine  the 
accuracy  of  fitting  shock  waves  with  a  full  Reynolds-averaged  Navier-Stokes  model  was  demonstrated  for  the  first 
time.  We  start  with  a  detailed  description  of  the  computational  procedure  developed  during  this  phase  of  our  work, 
followed  by  sample  calculations  and,  finally,  some  concluding  remarks  and  recommendations  for  future  work. 
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II.  DETAILS  OF  THE  COMPUTATIONAL  PROCEDURE 


1.  NON-DIMENSIONAUZATION 

Non-dimensional  physical  quantities  are  defined  as  follows: 

Let  p,p,a,f,S  be  dimensional  values  of  pressure,  density,  speed  of  sound,  temperature  and  entropy, 
respectively,  and  let  q  be  the  dimensional  velocity  vector.  Let  Cp  ,  Cy  be  the  specific  heats  at  constant  pressure 
and  constant  volume,  respectively,  and 

K  =  cr-cv,  y^cp/cy,  6  —  (y— 1)/2,  cp/H  =  y/26  (1.1) 

Assuming  the  values  of  pressure,  density  and  temperature  at  infinity  as  reference  values,  we  define  non- 
dimensional  pressures,  densities  and  temperatures  p,  p  and  T  as 

P  =  P/Poo,  p-P/pao)  T  =  f/Too  (1-2) 

Since 

P  —  'R-pT ,  Poo  =  TlpooToo  (13) 

the  basic  relation  between  p,  p  and  T  is  written  as: 

P  —  pT  (1.4) 

Also,  TlToo  is  assumed  as  the  square  of  the  reference  velocity,  qK fj  therefore,  the  non-dimensional  speed  of 
sound,  a,  is 

a  =  {yT)1'*  (1.5) 

A  suitable  reference  length,  zKf,  having  been  assumed,  the  reference  time,  i„{,  is  assumed  as  the  ratio  of 
the  reference  length  and  the  reference  velocity: 

f’rrf  =  xnt/qn{  (1.6) 

Finally,  we  assume  the  entropy  at  irfinity  equal  to  zero,  and  we  non-dimensionalize  the  entropy  by  the 

formula: 

S  =  S/(yK)  (1.7) 

The  Mach  number  at  infinity,  Moo,  is 

Moo  =  9co/a  00  (1.8) 

Additional  definition,"  are  needed  for  viscous  flows.  Let  A  be  the  divergence  of  q: 

A  =  V.q  (1.9) 


s  the  stress  tensor  minus  the  diagonal  tenns  containing  pressure,  Q  the  heat-flux  vector  and  $  the  dissipation 
function.  The  stress  tensor  (a  rather  complicated  function  of  space  derivatives  of  the  velocity  components) 
depends  linearly  on  the  visccsity,  p: 

s  =  '2p[def  q  -  -jAD]  (1-10) 

0 
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where  D  is  the  unit  diagonal  tensor  and 


defq=~[Vq  +  (Vq)T] 


(Ml) 


(the  apex  T  denotes  the  transpose  of  the  tensor).  The  heat-flux  vector  is  proportional  to  the  gradient  of  the 
temperature: 

Q  =  -iVf  (1.12) 

The  coefficient,  k,  is  the  thermal  conductivity.  We  also  define  a  “dissipation”,  $,  as 


$  =  2£(def  q  •  def  q  -  -A2) 

O 


(1.13) 


We  choose  reference  values  for  ji  and  k,  for  example  jioo  and  £<»  and  define  a  Reynolds  number, 


9oo  Poo  z  ref 


(1.14) 


and  a  Prandtl  number, 


The  continuity  equation,  in  dimensional  form,  is: 


(1.15) 


^+pV-q  =  0 


(1.16) 


The  momentum  equation,  in  dimensional  form,  is: 


Dq  1  1  . 

— 4  +  -r Vp=  -V  -s 
Dt  p  p 


(1.17) 


and  the  Lagrangian  increment  in  entropy  is: 


=  -L  [$  _  V  •  Q] 

Dt  pT  1 


(1.18) 


The  above  definitions  allow  the  three  equations  of  motion  to  be  recast  in  non-dimensional  form.  The 
continuity  equation  remains  formally  the  same,  without  the  bars: 


Dp 

^+pV.q  =  0 


(1.19) 


The  momentum  equation  becomes: 


where 


V  -  i  V-s 

Vm  ”  p  V  8 
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(s  is  formally  the  same  as  i  without  the  bars).  The  entropy  equation  becomes 


(1.22) 


v-mw:\ (I'23) 

where  the  non-dimensional  dissipation  and  the  non-dimensional  heat-flux  are  formally  the  same  as  $  and 
Q,  without  the  bars.  We  will  consider  the  viscosity  and  the  heat  conduction  coefficients  as  functions  of  the 
temperature: 

H  =  T3'2-1-*  (1.24) 

T  +  To/Toc  v  ’ 

(where  To  is  an  appropriate  constant;  for  example  110®  Kelvin),  and 


k^T3'* 


(1.25) 


2.  UPWIND  (“A”)  REFORMULATION  OF  THE  EQUATIONS 
The  equations  of  motion  (1.19),  (1.20)  and  (1.22)  can  now  be  reformulated  to  emphasize  the  role  of  the 
propagation  of  signals  in  the  convective  terms  and  the  diffusive  role  of  the  viscous  terms.  The  new  basic 
variables  are  a,  S  and  q.  From  the  definition  of  entropy  in  a  perfect  gas, 


S  =  — — pnp-  7lnp]  =  -—In - In  p 

2 y  '  n  2f6  7  7  K 

it  follows  that,  for  any  derivative,  here  denoted  by  a  prime, 

S'  =  ±  xt--L\ 

27  <  p  7  p 

7  6  a  7  p 

s>  =  l?L_l£ 

6  a  7  p 


Using  (2.3),  (1.19)  becomes 


1  Da  DS 

+  a  V  •  q  =  7a-—— 
6  Dl  H  '  Dt 


Using  (2.4),  (1.4)  and  (1.5),  (1.20)  becomes: 

~  +  jVa  -  a2VS  =  Vm  (2.6) 

To  make  the  computational  technique  for  viscous  flows  compatible  with  the  A-formulation  for  inviscid  flows, 
we  split  S  in  (2.5)  into  two  terms: 

5  =  (i  +  ^)S  (2.7) 
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recasting  (2.5)  into  the  form: 


At  this  point,  (2.8),  (2.6)  and  (1.22)  can  be  treated  by  the  general  procedure  exposed  in  [1],  with  the  addition 
of  the  source  terms,  26a Vt,  Vm,  and  Vt,  respectively. 

The  flow  is  two-dimensional,  and  the  computational  grid  is  orthogonal.  In  the  initial  phase  of  the  computa¬ 
tion,  the  rigid  walls  and  the  computational  frame  move  from  right  to  left,  accelerating  from  a  state  of  rest 
to  the  constant  velocity, 

?oo  =  Moov/7  (2.9) 

.  The  equations  of  motion  are  recast  for  the  general  case  of  an  accelerating  frame.  This  changes  only  the 
momentum  equation  (2.6),  to  which  a  term,  ?oo*I  must  be  added  in  the  right-hand  side,  with  q^t  being  the 
acceleration  of  the  frame  and  I  the  unit  vector  of  the  x-axis. 

We  rewrite  the  equations  of  motion  after  splitting  the  Lagrangian  derivatives  into  local  derivatives  and 
convective  terms: 

1%  +  q  ‘  T  +  aV  ‘  q~  “(§T  +  q  -  vs)  =  26aV>  (210) 

^  +  (q  •  V)q  +  |Va  -  a2V5  =  qxtl  +  Vm  (2.11) 

f  +  q-VS  =  K  (2.12) 

These  equations  are  the  same  as  equations  (1.2)  in  [1],  with  the  addition  of  the  source  terms  in  the  right-hand 
sides. 

We  proceed  as  explained  in  Section  6  of  (1],  Four  equations  are  obtained  containing  two  unit  vectors, 
orthogonal  to  each  other,  i  and  j,  that  we  choose  parallel  to  the  coordinate  lines  at  each  point  (these  vectors 
are  a  particular  case  of  the  unit  vectors  n  and  r  of  [1]).  Using  indices  to  denote  partial  derivatives  and 
calling  u  and  v  the  components  of  q  along  the  coordinate  lines,  the  equations  (similar  to  (2.1)  of  [1],  but 
with  right-hand  sides),  are: 

at/6  +  -  aSt  +  (q  +  ai)  •  [V(a/5  +  u)  -  aVS]  -f  aj  •  Vv  -  0v  +  F  =  2 6aV,  +  ?oo<I  •  i  +  Vm  •  i 

at/6  -  «(  -  aSt  +  (q  -  ai)  •  [V(a/£  -  u)  -  aVS]  +  a  j  -Vv  +  fiv  +  F  =  2  6aV,  -  gTO(  I  •  i  -  Vm  •  i 

at/6  +  vt  -  aSt  +  (q  +  aj)  •  [V(a/6  +  v)  -  aVS]  +  ai  •  Vti  + /?u  -f  F  =  26a V,  +  q^l  ■  j  +  Vm  ■  j 

at/6  -  v(  -  aSt  +  (q-  aj)  •  [V(a/<5  -  v)  -  aVS]  +  ai  •  V«  -  /?u  +  F  =  2 6aV,  -  q^,!  •  j  -  Vm  •  j 

(2.13) 

with 

a  =  arccos(I  •  i),  k  =  ixj  (2.14) 

P  -  qVa,  F  =  akxqVa  (2.15) 

Adding  the  four  equations  (2.13),  and  subtracting  (2.10)  multiplied  by  2,  we  obtain: 

(2 /6)at  -  2 aSt  +  2aq  ■  VS  +  Ux  +  U2  +  U3  +  UA  -  (2/6)q  ■  Va  =  46aVt  (2.16) 
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Subtracting  the  second  of  (2.13)  from  the  first,  we  obtain: 

2u, +  i/i  -  f/2  =  2g00<cosa+ 2Vm -i  (2-17) 

Subtracting  the  fourth  of  (2.13)  from  the  third,  we  obtain: 

2v(  +  Ua  -  U4  =  2goo*  sin  a  +  2Vm  -j  (2.18) 

where 

Ui  =  (q  +  ai)  •  [V(a/£  +  u)  —  vVa  -  aVS] 

U2  =  (q  —  ui)  •  [V(a/6  —  u)  +  vVa  -  aVS] 

Ua  =  (q  +  aj)  •  [V(a/5  + 1>)  +  uVar  -  aVS]  (2-19) 

U4  =  (q  -  aj)  •  [V(a/6  -  r)  -  uVa  -  a  VS] 

These  equations  are  the  counterparts  of  (2.2),  (2.3)  and  (2.4)  in  [1],  with  a  different  grouping  of  terms  to 
eliminate  the  explicit  appearance  of  F  and  /?.  Note  that  the  first  parenthesis  in  the  right-hand  sides  of  (2.19) 
contain  vectors.  For  example,  in  the  first  equation, 

q  +  ai  =  (u  +  a)i  +  vj 

When  expanding  the  expressions,  the  terms  affected  by  u  +  a  must  be  approximated  using  differences  taken 
from  the  side  originating  the  (u  +  a)-signal.  Such  terms  cannot  be  simplified  with  similar  terms  in  the 
second  equation,  affected  by  the  coefficient  u  -  a.  Contrariwise,  the  terms  affected  by  the  coefficient  v 
can  be  simplified  between  the  first  two  equations,  because  both  have  the  same  domain  of  dependence.  The 
important  simplifications  that  the  above  argument  implies  are  shown  in  what  follows,  in  the  case  of  orthogonal 
coordinates. 

The  physical  plane  is  defined  by  a  complex  variable  2  =  x+iy  and  is  conformally  mapped  onto  another  plane, 
defined  by  a  complex  variable  (  =  (  +  »'»?•  The  grid  in  the  (-plane  is  Cartesian;  therefore  its  counterpart 
in  the  2-plane  is  orthogonal.  It  may  be  necessary  to  stretch  the  coordinate  lines  in  the  (-plane  to  obtain  a 
better  resolution  in  the  2-plane.  In  this  case,  the  rectangular  grid  in  the  (-plane  is  non-conformally  related 
to  another  rectangular  grid  in  a  plane,  defined  by  Cartesian  coordinates  X  and  Y,  where  X  is  a  function  of 
(  alone  and  Y  is  a  function  of  77  alone;  this  does  not  impair  the  orthogonality  of  the  final  mesh.  All  gradients 
appearing  in  the  left-hand  sides  of  (2.16),  (2.17)  and  (2.18)  are  conveniently  computed  in  the  (-plane.  With 

9  =  ~j~i  G  =  M.  <f>  =  <f>\  +  i<t>2  =  ~  (2.20) 

we  have: 

Va  =  G(afi  +  a^j)  =  —G(<j>  2i  +  <^jj)  (2-21) 

and 

U\  =  G(u  +  a)[(aj6  +  u)^  -  -  aS^]  +  Gv\(a/f>  +  u)n  -  -  aS,] 

U2  =  G(u  -  a)[(a/o  -  u)(  -f  va(  -  aS^]  +  Gv[(a/6  -  u)n  -f  van  -  aS,,] 

f/3  =  Gu[(a/6  +  v){  +  uoj  -  aS{]  +  G(r  +  a)[(a/6  +  v),  +  -  a5,,} 

b\  =  Gu[(a/6  -  v){  -  uct{  -  aS^]  -f  G(v  -  a)[(a/6  -  u),,  -  ua^  -  a.Sp] 
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q  •  [V(a/S)  -  aVS]  =  G[(ua(/6  -  auS()  +  (t >av/6  -  avSn)] 


(2.23) 


The  above  mentioned  simplifications  occur  in  (2.16)  because  the  second  terms  in  U\  and  f/a  and  the  first 
terms  in  Uz  and  U4  cancel  out  the  right-hand  side  of  (2.23).  On  the  other  hand,  the  second  terms  of  U\  and 
C/2  reduce  to  2Gv(un  -  van),  and  the  first  terms  of  Uz  and  U4  reduce  to  2Gu(v{  +  uazi)  when  computing 
U\  —  Uz  and  Uz  —  U4,  respectively.  For  better  clarity  in  coding,  we  define 

Af  =  G(u  +  a),  A f=G(u-a),  Af  =  Gu 


Af  =  G(v  +  a),  =  G(v  -  a),  Aj  =  Gv 

R*  =  a/6  +  u,  R$  =  a/6  -  u,  R%  =  a/6  -ft),  R%  =  a/6  —  v 


(2.24) 


ft  =  -A i(Ri(  -  va(  -  aS(),  f*  =  -Af  (R^  +  va(  -  aSf ),  / f  =  -2Af  (vf  +  ua€) 
fi  =  -Ax' (R\„  +  uav  -  aSn),  =  -A2  (/^  -  tiarn  -  aS,,),  /3y  =  ~2A£ (u,  -  va„) 


(2.25) 


The  equations  of  motion  are  now: 


/f*-A?s«,  /r=-Ars„ 


s<  =  /f-f/y  +  v. 


at  —  2^iX  +  /*  +  /y  +  /I"  +  2aSt  -f  460!^] 
=  ^[fi  ~  f  *  +fz+  2?oo«  cos  a  -f  2Vm  •  i] 
*>t  =  ^[/iy  ~  fz  +  fz  ~  2?oot  sin  a  +  2Vm  •  j] 


(2.26) 


(2.28) 


(2.29) 


(2.30) 


We  must  now  evaluate  the  viscous  terms,  V,  and  Vm.  In  orthogonal  coordinates,  generated  by  conformal 


mapping, 


en -=  G(ti{ -f  vfo),  «22  =  G{vn  -  u<f>i),  ei2  =  —  (v{  +  u,  -  ufc  -f  tx£i) 


A  =  en  -f  e22 


The  two  components  of  V  •  s  are 


(2.31) 

(2.32) 


2G{[fi(el\  -  - A)](  +  (pe !2)n  +  2pei2<£2  +  /i(e22  -  ejj)^>i]} 


(2.33) 


2G{(fie\z)(  -f  (p(e22  -  -A)],  -  2pei2</>)  -f/i(en  -  e22)^2]} 


Using  a  classical  notation,  we  define: 


(2.34) 


T11  =  2/i(en  -  ^A),  r12  r=  2/iei2,  r22  =  2p(e22  -  ^A) 
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(2.36) 


I 

I 

I 


I 

I 

I 


Therefore,  the  two  components  of  Vm  can  be  written  as: 

vm  •  i  =  ”v^—  [(m)e  +  (Ti2)>;  +  2n2^2  +  (^22  -  m)&] 

rCe  p 

vm  •  j  =  -  —  [(n2)e  +  (722)11  “  2ri2^i  +  (rn  -  r22)^2] 

p 

From  its  definition  (1.13), using  (2.31)  and  (2.32),  $  can  be  written  as: 

*  =  4/»[^(cii  +  ej2  -  ene22)  +  c?2] 

Finally,  the  term  V  •  Q  is  evaluated  as 

V  •  Q  =  -V  •  (kVT)  =  -G2[(kT()(  +  (kTn)v) 

Let  uc  and  vc  be  the  Cartesian  components  of  q;  thus, 

u  =  uc  cos  a  +  vc  sin  a 
v  =  —uc  sin  a  +  vc  cos  a 

and 

uc  =  u  cos  or  —  v  sin  a 
vc  =  u  sin  a  +  v  cos  o 


(2.37) 

(2.38) 

(2.39) 

(2.40) 

(2.41) 


A)  Left  computational  boundary  -  The  flow  is  generally  considered  inviscid  to  the  left  of  the  computa¬ 
tional  boundary.  The  observer  moves  with  the  rigid  walls;  at  the  end  of  the  accelerating  phase,  the  flow  at 
infinity  appears  as  moving  towards  the  observer  at  the  speed  defined  by  (2.9).  If  Moo  is  greater  than  1,  a 
is  forced  equal  to  ^/y,  uc  to  g^, ,  vc  and  S  to  zero  over  the  entire  left  computational  boundary  after  each 
updating  defined  by  (2.27),  (2.28),  (2.29)  and  (2.30),  and  the  values  of  /*,  fj  ( i  =  1,2, 3,4)  are  all  defined 
using  (2.40)  before  such  updatings. 

If  Mto  is  less  than  1,  vc  is  still  assumed  to  vanish  along  the  left  computational  boundary  (this  is 
tantamount  to  directing  the  flow  horizontally  by  means  of  an  infinite  cascade  of  vanes,  and  it  is  not  a 
disturbing  assumption).  The  values  of  all  /,y ’s  we  defined  accordingly.  The  steady  flow  conditions  to  the 
left  of  the  boundary  are  expressed  by  the  Constance  of  the  stagnation  speed  of  sound  and  by  the  vanishing 
of  the  entropy.  Such  conditions  are  sufficient  to  determine  the  values  of  a  and  u,  separately,  when  combined 
with  the  value  of  R$ ,  that  is  obtained  correctly  using  inside  information  only.  Specifically, 

do  =  a2  +  S(u2  +  v2)  =  a2  -f  6u2/  cos2  a  (2-42) 

and 

1$  =  a/S  -  v  (2.43) 

From  (2.42)  we  obtain: 

aat/5  +  uu(/ cos2  a  =  0  (2-44) 
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and  using  (2.28)  and  (2.29)  without  viscous,  acceleration  and  entropy  terms,  f*  can  be  defined  as  follows: 

fx  =  [(u/  cos2  a  -  a)f*  -  ufz  /  cos2  a  -  a{f]'  +  f\ )]/(«/  cos2  a  +  a)  (2.45) 

In  particular,  if  the  left  boundary  is  vertical, 


fi  = 


ti  +  a 


(2.46) 


Resetting-  The  above  procedure  is  analytically  correct;  it  may  produce  disturbing  effects  in  a  numerical 
computation,  if  continued  for  thousands  of  steps.  In  fact,  (2.45)  provides  a  correct  evaluation  of  the  time 
derivatives  of  a  and  u  at  the  left  boundary  (that  is,  an  evaluation  assuming  that  a0  is  constant  in  time  and 
space).  The  updating  of  a  and  u  at  the  left  boundary  according  to  (2.45),  however,  is  affected  by  truncation 
errors  and,  in  a  long  run,  such  errors  may  accumulate  resulting  in  minor,  but  not  insignificant,  variations 
in  the  effective  value  of  ao-  The  inconvenience  can  be  overcome  by  resetting  a0  to  its  exact  value  after  each 
updating,  and  recomputing  a,  u  and  v  accordingly.  From  (2.42)  and  (2.43),  ti  turns  out  to  be  defined  by  the 
equation: 

Au*  +  2Bu  +  C  -  0  (2.47) 

with 

A  —  6(6  +  1/cos2  a),  B  =  62Rf,  C  =  6\R*)2  -  a2 
[the  +  sign  must  be  used  in  front  of  the  square  root  in  solving  (2.47)];  then,  a  follows  from  (2.43)  and 


v  =  —  u  tan  a  (2.48) 

.  When  the  flow  is  supersonic  or  resetting  is  applied,  the  derivatives  computed  on  the  left  boundary  are  not 
used  to  integrate  the  equations  of  motion  on  the  left  boundary,  but  the  values  of  ff's  and  f^'s  must  be 
evaluated  at  the  first  level  of  a  two-level  scheme,  because  they  are  needed  at  the  second  level,  as  seen  later 
on. 

B)  Upper  boundary  -  We  distinguish  two  geometrical  possibilities.  If  the  upper  computational  is 
horizontal,  it  will  be  considered  as  a  rigid,  inviscid  wall,  along  which  =  f\  =  0  and  is  given 
by 

ft  =  ft  (2-49) 

to  satisfy  the  inviscid  boundary  condition,  v  =  0. 

If  the  upper  computational  boundary  is  curved,  with  the  concavity  upwards,  it  is  considered  as  an  entry 
boundary  for  the  flow.  Therefore,  if  the  flow  is  supersonic,  fill  values  arc  prescribed  and  kept  invariant  in 
time  at  the  end  of  the  accelerating  phase.  If  the  flow  is  subsonic,  a  routine  similar  to  the  one  explained  for 
the  subsonic  left  boundary  will  be  applied,  interchanging  the  roles  of  u  and  v.  Specifically, 

fi  =  [(v  +  o  sin2  o)fl  +  vf£  +  a(/j*  +  /* )  sin2  a]/(t>  -  a  sin2  a)  (2.50) 

and,  for  resetting,  (2.47)  can  be  used,  with 

.4  =  6(5  sin2  a  -f  1),  B  =  -67R]  sin2  a,  C  =  [<52(tff  )2  -  ajjjsin2  a 
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C)  Lower  boundary  -  The  lower  boundary  is  composed  of  two  parts:  a  free-stream  region  to  the  left  of 
the  leading  edge  of  the  wall,  and  the  wall  itself.  The  former  is  always  a  horizontal  segment  and  it  is  consider 
fully  inviscid.  Therefore,  it  acts  as  an  inviscid  rigid  wail,  with  t;  vanishing  identically,  and  the  only  numerical 
condition  is 

/r = i\  (2.5i) 

At  any  point  on  the  wall,  u  =  v  =  0.  If  the  wall  is  isothermal,  its  temperature,  Tw  is  prescribed.  If  it  is 
adiabatic,  Tw  equals  the  computed  value  of  T  at  the  grid  point  just  above.  The  pressure  is  also  considered 
equal  on  any  wall  point  and  the  point  above  it,  in  any  case.  Therefore,  if  the  wall  is  isothermal,  a  follows 
from  its  definition  (1.5)  and  5  from  the  Constance  of  p,  using  (2.4)  and  (1.5): 

Sr,  -  an/{Sa)  (2.52) 

If  the  wall  is  adiabatic,  Sn  =  a,  s=  0.  Thus,  all  quantities  at  the  wall  are  easily  known  and  no  computation 
of  //*’ s  or  viscous  terms  is  necessary.  In  a  two-level  technique,  ,f%  and  f\  at  the  wall  are  needed  at 
the  second  level  for  the  grid  row  immediately  above  the  wall  to  provide  a  second-order  correction,  since  v  is 
normally  positive.  They  may,  however,  be  taken  equal  to  their  values  on  the  row  above  the  wall. 

D)  Right  boundary  -  On  the  right  boundary,  a  zero-th  order  extrapolation  for  all  variables  can  be  used. 
In  a  supersonic  zone,  the  extrapolation  is  obviously  harmless,  because  no  signal  propagates  upstream.  Inside 
the  subsonic  zone  of  a  boundary  layer,  the  approximation  is  acceptable  since  the  flow  is  essentially  parabolic. 
The  only  dubious  case  occurs  with  a  practically  inviscid  subsonic  boundary,  where  f*  at  points  immediately 
next  to  the  boundary  is  affected  by  the  values  of  a  and  u  at  the  boundary  itself.  In  the  computation,  however, 
the  error  can  be  minimized  by  letting  f*  =  0  both  at  the  boundary  and  on  the  grid  column  immediately 
before  it. 

3.  OUTLINE  OF  THE  COMPUTATIONAL  CODE,  IN  THE  ABSENCE  OF  SHOCKS 

In  the  absence  of  shock,  the  computational  code  is  structured  as  follows: 

The  main  program  calls  a  subroutine  ELDATA  that  reads  the  input  data,  necessary  to  define  the  geometry 
of  the  rigid  walls  and  to  generate  a  proper  grid,  the  value  of  Moo,  and  the  Reynolds  and  Prandtl  numbers 
according  to  their  definitions  (1.14)  and  (1.15).  Other  data  are  parameters  to  define  the  coordinate  stretchings 
between  £  and  X  and  between  rj  and  Y. 

Then  it  calls  a  subroutine  GRID  that  generates  the  grid  and  provides  double  arrays  with  the  values  of 
G,  <f> l,  4> 2,  cos cr  and  sina,  and  single  arrays  with  the  values  of  X(  and  Yv.  The  values  of  G  are  obtained 
numerically  as 

G  =  jA(/Az|  (3.1) 

using  centered  differences.  Similarly,  a  is  obtained  as 

q  =  arctan(OAz/9?Az)  (3.2) 

using  centered  differences.  Finally,  <f>,  and  4>2  are  obtained  as 

<h=G</G,  fo  =  -Gn/G  (3.3) 
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and  the  derivatives  of  G  are  approximated  numerically,  using  centered  differences,  evaluated  on  the  ( X ,  Y)~ 
plane  and  interpreted  in  terms  of  t  and  17  via  X(  and  Vi,. 

The  field  is  initialized  as  a  gas  at  rest,  with  non-dimensionalized  values  of  the  thermodynamical  pa» 
rameters  as  mentioned  above  for  the  conditions  at  infinity  ( T  —  1,  a  =  ftj,  5  =  0).  Then  the  sequence  of 
computational  steps  is  started. 

In  it,  qoot  is  computed;  a  subroutine  VISC  is  called  to  compute  V,  and  Vm;  next,  (2.25),  (2.26),  (2.27) 
and  (2.28)  are  integrated  at  the  predictor  level  in  a  subroutine  PRED.  The  entire  procedure  is  repeated, 
but  PRED  is  replaced  by  a  subroutine  CORR  for  the  integration  at  the  corrector  level.  This  completes  one 
computational  step.  Boundary  conditions  are  enforced  by  calling  a  subroutine  BOUID  both  from  PRED  and 
from  CORR. 

We  outline  now  the  four  subroutines  used  in  each  step. 

PRED  contains  a  loop  to  define  A* .  Af  («  =1,2, 3),  and  R? ,  Rf  (i  =  1,2).  The  time  step  size  is  evaluated 
according  to  the  CFL  rule;  if  necessary,  a  stability  coefficient,  less  than  1,  is  used  to  reduce  it.  A  second 
loop  computes  ft ,  ft (i  =  1,2, 3, 4, 5)  using  averaged  values  of  A-s  over  two  adjacent  points  and  two-point 
differences,  both  from  the  side  from  which  the  pertinent  signal  comes.  In  defining  the  /- s,  each  one  is  divided 
by  8.  A  factor  of  2  is  due  to  averaging  the  A-s;  another  factor  of  2  is  the  one  in  front  of  the  brackets  in 
(2.26),  (2.27)  and  (2.28);  the  third  factor  of  2  makes  the  updating  take  place  over  a  half  time  step  only. 

Successively,  the  boundary  conditions  are  enforced.  Then,  a  third  loop  defines  St,  at,  u,  and  vt,  in  that 
order,  using  the  formulas: 

St  =  -51/,  +  /*  +  ft 

a<  =  &(Ji  +  ft  +  fx  +  ft  +  aSt  +  6al/,) 

=  ft  ~  ft  +  ft  +  -5(9001  cos  a  +  Vu)  (3.4) 

=  ft  -  ft  +  ft  ~  -5(9oo«  sin  a  —  Vv) 

and  updates  5,  a,  u  and  v.  Resetting  of  all  boundary  values  is  made.  In  a  final  loop,  all  the  f-s  are  stored 
in  “old”  arrays,  /^,  etc.  The  ft  of  the  left  boundary  are  also  stored  in  the  “old”  arrays  for  the  extra  row 
to  the  left;  similarly,  the  ft  of  the  upper  boundary  are  stored  in  the  “old”  arrays  for  the  extra  row  above 

it.  Such  values  will  be  used  in  CORR  to  cor:*  *.t  values  on  the  boundaries  and  are  acceptable  in  all  supersonic 

cases  on  the  left  boundary.  On  the  upper  boundary,  if  curved,  Af  is  generally  positive;  therefore,  ft  must 
be  corrected  in  BOUID.  Similarly,  / *  must  be  corrected  in  BOUID  on  the  left  boundary,  if  the  flow  is  subsonic 

CORR  follows  the  general  pattern  of  PRED.  The  time  step,  however,  is  no  longer  evaluated.  All  the  f-s 
are  now  divided  by  4  only,  and  their  “old"  counterparts,  computed  in  PRED,  are  subtracted  from  them.  Note 
that  such  “old”  values  must  be  taken  at  the  neighboring  point,  from  the  side  fro  which  the  signal  proceeds 
No  storage  of  f-s  into  “old”  arrays  is  now  necessary. 

VISC  begins  by  defining  T  as  a  function  of  a  at  all  grid  points,  and  by  linearly  extrapolating  u,  v,  a,  T 
and  S  at  the  grid  column  outside  the  computational  field,  to  the  right  of  the  right  boundary  This  is  made 
to  enable  centered  differences  to  be  taken  at  boundary  points,  using  the  same  scheme  as  at  interior  points 
For  the  same  reason,  u,  a  and  T  are  defined  symmetrically,  and  t  is  defined  antisymmetrically  at  the  lower 
boundary  in  front  of  the  wall.  Along  the  wall,  u  and  v  are  defined  antisymmetrically,  whereas  a,  T  and  5 
are  extrapolated  linearly. 
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Then,  at  all  points  except  at  the  upper  boundary,  centered  differences  are  used  to  approximate  ti(,  v^, 
T(,  u,,  v,  and  Tn\ ft  is  computed  as  a  function  of  T  according  to  (1.24)  and  (2.31),  (2.32),  (2.35)  and  (2.38) 
are  evaluated.  Using  a  function  k(T)  as  defined  in  (1.25),  Qi  =  GkJ\  and  Qj  =  GkT^  are  computed.  Again, 
Tn,  ri2  and  tJi  are  linearly  extrapolated  outside  the  right  boundary,  t22  and  Q2  aie  made  symmetric,  and 
is  made  antisymmetric  across  the  wall. 

Finally,  at  all  points  except  on  the  lower  and  upper  boundaries  centered  differences  are  used  to  evaluate 
(tu)(,  (ru){,  (*?i){.  (Tujq,  (rJ2 )fj  and  (Qj)n.  After  this,  all  the  elements  needed  to  build  V,  and  Vm  are 
available.  The  viscous  terms  are  made  equal  to  zero  at  all  grid  points  to  the  left  of  the  grid  line  passing 
through  the  leading  edge  of  the  wall. 

BOUID  redefines  the  /’ s  at  the  boundaries  when  necessary.  As  mentioned  before,  no  correction  is  nec¬ 
essary  on  the  left  boundary,  if  the  flow  is  supersonic.  If  it  is  subsonic,  f*  is  redefined  according  to  (2.45). 
On  the  upper  boundary,  is  always  redefined  according  to  (2.50),  if  the  boundary  is  curved.  Otherwise, 
(2.49)  is  used.  At  the  right  boundary,  /*  is  set  equal  to  zero  if  the  flow  :i  subsonic.  At  the  lower  boundary, 
to  the  left  of  the  wall,  (2.51)  is  used.  On  the  wall,  fi  =  fi  ~  Vm  •  j/2;  this  term  is  only  used  to  produce  a 
value  of  to  be  used  in  CORK  at  the  grid  row  above  the  wall. 

4.  COMPUTATION  OF  SHOCKS 

Shocks  arc  considered  as  sharp  discontinuities.  Each  shock  is,  thus,  a  line  separating  a  low-pressure  region 
from  a  high-pressure  region.  At  each  point  on  a  shock,  a  normal  to  the  shock  can  be  characterized  by  a  unit 
vector,  N.  All  values  on  the  low-pressure  side  of  a  shock  will  be  denoted  by  an  index,  A,  all  values  on  the 
high-pressure  side  by  an  index,  B. 

In  what  follows,  it  is  convenient  to  consider  three  set  of  unit  vectors  at  each  shock  point: 

1)  I  and  J,  parallel  to  the  Cartesian  axes  x  and  y, 

2)  i  and  j,  parallel  to  the  grid  lines  (  and  rj, 

3)  N  and  T,  normal  and  parallel  to  the  shock. 

According  to  (2.14), 

I  ■  i  =  coe  a,  I-j  =  -sina,  Ji  =  sina,  Jj  =  cosa  (4.1) 

Let  V’  be  the  angle  between  N  and  I,  8  the  angle  between  N  and  i.  Then, 

0  =  V'-Q  (4.2) 

The  two  components  of  N  along  i  and  j  will  be  called  N\  and  N 2,  respectively: 

.V]  =  N  •  i  =  cos  8,  N2  =  N  •  j  =  sin  8  (4.3) 

The  velocity  vector,  q,  can  be  decomposed  along  the  three  frames  defined  above: 

q  =  ucI  +  rcJ  =  ui  +  uj  =  uN  -f  fT  (4.4) 

Consequently, 

ii  =  tiWi+  vA'j,  vzz-uN2  +  vN\  (4.5) 
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and,  conversely, 


(4.6) 


u  =  uNi  —  vNt,  v  =  uNj  +  vNi 


The  motion  of  a  shock  is  defined  by  the  velocity  (W)  of  each  of  its  points  along  the  local  normal.  The 
expression: 


(4.7) 


is  the  normal,  relative  Mach  number  of  the  shock  at  each  point;  here  it  will  be  called  the  shock  Mach  number, 
for  brevity.  The  Rankine-Hugoniot  conditions  can  be  written  in  the  form: 


[(■yA/2  -  £)(1  +  A  Af2)]1/2 
as  -aA  (1  +  6)M 


+  aA 


VB  =  V.\ 


1  -  Af2 

(1  +  6)M 


„  „  .  1  ,  7M2-6  ,  (l+6)M7 

SB=s^2h  bnr+r“7toT+M#r 


(4.8) 


We  define 


E  =  (as  +  6\ us  -  uA\)/aA 


(4.9) 


From  these  equations  we  obtain 

E  =  (1  -M)M  [v^(7A/2  ~  ^)(1  +  6M^  +  ^  (410) 

which  yields  the  basic  relationship  between  E  and  M ;  consequently,  if  E  is  known,  M  can  be  computed  from 
(4.10).  For  computational  purposes,  it  is  convenient  to  use  a  relationship  between  an  increment  in  M  and 
the  corresponding  increment  in  E: 


AE  = 


AM  2j6M*  +  y-67 

l  +  «  ^A^-axi  +  ^M2) 


o+o§] 


(4.11) 


In  the  framework  of  the  computational  technique  exposed  in  Section  2,  we  observe: 

a)  once  we  have  the  basic  variables,  a,  5,  u  and  5  on  the  low-pressure  side  of  the  shock,  and  the  shock  Mach 
number,  the  shock  problem  is  solved;  in  fact,  (4.8)  provide  the  values  on  the  high-pressure  side  and  (4.7) 
provides  the  velocity  of  the  shock; 

b)  the  values  of  a,  S,  u  and  v  on  the  low-pressure  side  of  the  shock  depend  only  on  upwind  values  and  are 
accurately  defined  using  the  A-scheme  with  upwind  differences  only; 

c)  knowledge  of  the  direction  of  the  normal  to  the  shock  is  imperative  to  get  u  and  C  from  u  and  t>,  and  to 
get  u  and  v  on  the  high-pressure  side  of  the  shock  from  ti  and  5;  moreover,  it  is  necessary  to  move  the  shock 
point  correctly,  once  W  has  been  determined; 

d)  the  shock  Mach  number  is  obtained  from  E,  using  (4.10)  or  (4.11),  but  E  must  be  obtained  from  the  low- 
and  high-pressure  side  of  the  shock. 
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By  writing  (4.9)  in  the  form: 


(4.12) 


E  =  6[aB/6  -  iifl  +  uA}/aA 

and  comparing  with  the  last  of  (2.23),  we  see  that  (4.12)  can  be  written  as: 

X  =  6(R7B  +  uA)/aA  (4.13) 


where  R7B  is  to  be  computed  using  downstream  information,  since  the  flow  in  the  high-pre3sure  side  is 
subsonic,  relatively  to  the  shock.  The  evaluation  of  E  is  thus  feasible  without  using  differences  across  the 
shock. 

Clearly,  the  most  difficult  step  in  fitting  a  shock  consists  of  determining  the  direction  of  the  normal. 
When  a  shock  point  is  detected,  we  consider  the  jump  in  q  between  the  two  points  bracketing  it,  and  define 
the  normal  by  its  unit  vector  as 


N  = 


A6  fq 
\Abfq\ 


(4.14) 


The  same  procedure  has  not  been  proven  generally  safe  to  evaluate  the  normal  at  a  shock  point  that 
is  part  of  an  existing  link  of  shock  points.  In  this  case,  we  prefer  to  use  a  more  cumbersome  procedure. 
For  a  given  shock  point,  located  on  the  m-th  horizontal  row  to  the  right  of  the  n-th  vertical  column,  we 
explore  some  intervals  to  the  left  and  the  right  of  it  on  both  the  m  -  1  and  the  m  -f  1  row.  If  no  points 
are  found,  the  shock  point  is  considered  isolated  and  it  is  dropped.  If  some  point  is  found,  the  slope  of  the 
shock  is  computed  as  the  slope  of  the  line  joining  the  geometric  center  of  the  points  on  the  lower  row  and 
the  geometric  center  of  the  points  on  the  upper  row,  if  the  flow  is  subsonic  downstream  of  the  shock.  If  the 
shock  is  supersonic  downstream,  the  line  is  drawn  between  the  shock  point  itself  and  the  geometric  center 
of  the  points  on  the  row  that  contains  the  “upstream”  shock  point. 

The  shock  is  computed  only  at  the  end  of  each  computational  step,  after  all  the  values  at  grid  nodes 
have  been  updated  by  predictor  and  corrector.  The  values  at  points  A  and  B  on  the  shock  itself  cannot 
be  updated;  it  is  indeed  in  the  nature  of  a  numerical  scheme  on  a  discreet  grid  to  let  all  points  that  are 
not  nodes  remain  undefined.  A  simple  way  to  circumvent  the  difficulty  is  to  assume  that  such  values  are 
the  same  as  the  values  at  the  nearest  node  on  either  side  (0-th  order  extrapolation).  The  approximation 
is  generally  sufficient  for  unsteady  flows;  it  may  prevent  convergence  of  residuals  to  machine-zero  in  steady 
flows  asymptotically  reached,  if  a  shock  point  happens  to  lie  too  close  to  a  grid  node.  The  shock  point  may 
periodically  move  from  one  grid  cell  to  the  next,  and  back  again;  at  each  transition,  the  extrapolated  values 
may  be  sufficiently  different  to  create  a  local  disturbance  that  propagates  downstream  in  successive  steps. 
The  difficulty  is  eliminated  by  using  a  more  sophisticated  extrapolation  technique  to  define  point  A  and 
point  B.  For  example,  for  point  A,  located  to  the  right  of  grid  p^int  n,  we  can  define: 


i  —  |-j  Zn|/|ln+l  —  In  | 

^  1  —  —  Zn |/]lii  ~  In—  1 1  (4.15) 

<2  =  |->  ~  Zn-l|/|ln-l  -  Zn-2| 

and  then 

«1  =  On  +  <l{c>n  -  On-l) 

°2  =  an-l  +  <2(Qn_l  -  Un-2?  (4.16) 

aA  =  <aj  +  (1  -  c)a2 
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and  similarly  for  5,  u  and  v.  It  is  easy  to  see  that  the  values  at  A  do  not  change  if  the  shock  point  mo.'es 
from  the  right  boundary  of  cell  (n,n  +  1),  with  <  close  to  1,  to  the  left  boundary  of  cell  (n  +  l,n  +  2),  with 
<  close  to  0. 

The  values  at  A  must  now  be  considered  as  the  final  updated  values  for  the  computational  step  on 
hand.  The  values  at  B  still  have  to  be  updated  after  the  shock  itself  has  been  upaated.  Only  the  combined 
quantity,  Ri  is  acceptable  so  far.  We  proceed,  thus,  to  get 

Au  =  ( ub  -  u a)Ni  +  (vb  -  va)N3  (4.17) 

and 

E  =  (aB -M|Au|)/a„  (4.18) 

Now,  (4.11)  may  be  applied  to  get  the  increment  of  M  through  the  time-step  from  the  increment  of  E;  having 
M,  (4.8)  are  applied  to  get  new  values  at  B;  ub  and  ti/j  will  be  decoded  from  ub  and  vg  using  (4.6).  Final 
values  at  the  node  next  to  B  (n  +  1,  say)  are  obtained  from  formulas  such  as 


«n+i  =  as  +  *3(an+2  -  ob)  (4-19) 

with 

<3  =  |*n+l  -  2.  l/kn+2  “  2. 1  (4.20) 

The  shock  velocity  is  computed  from  (4.7).  The  displacement  of  the  shock  point  along  a  (-line  is  obtained 
as 

A(  =  W-^-At  (4.21) 

N  i 

la  general,  the  formation  of  a  shock  by  coalescence  of  characteristics  is  detected  as  follows: 

On  each  (-line,  the  first  difference  of  a  and  the  second  difference  of  A2  (briefly,  A)  are  approximated  at 
each  point  n  by 

Aa  =  on+j  -  a„_i  +  2(an+j  -  an  -  2)  (4.22) 

and 

A2A  =  ~A„+1  -  A„_2  -  An—  1  —  An  +  2(A„^  2  f-  An— 3)  (4-23) 

respectively.  These  expressions  have  been  obtained  by  replacing  the  data  with  their  linear  mean-square 
approximations,  to  eliminate  the  influence  of  minor  numerical  irregularities.  If  the  flow  goes  from  supersonic 
to  subsonic  (th’t  is,  if  A„  >  0  and  A„+l  <  0)  and  -An+i/An  >  1.1,  the  point  on  the  A'-line  where  the 
interpolated  A  vanishes  is  marked.  Other  points  marked  are  those  where  A2An  <  0,  A2An+1  >  0  and 
A  a  >  0.  If  any  of  such  points  is  found  along  the  A-line,  a  tentative  value  of  E  is  defined  using  <7„+i  and  q„ 
instead  of  ug  and  «m,  and  an+i,  an  instead  of  ag ,  a  a  The  point  is  assumed  as  a  new  sho'k  point,  provided 
that  E  >  1.03  and  no  other  shock  points  exist  on  the  cells  between  n  -  4  and  n. 

If  a  shock  is  genctated  by  a  wedge  or  a  sharp  corner  in  a  rigid  wall,  a  third-degree  equation  is  solved 
iteratively  to  toet  the  direction  of  the  normal  to  the  shock  at  its  point  of  origin  With  <7  being  tht  „quare  of 
the  sine  of  the  flow  deflection,  r  the  square  of  the  cosine  of  the  angle  between  the  normal  and  the  upstream 
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1 

I 


direction  of  the  flow,  Ma  —  va/i>a,  di  =  1/AfJ,  b  =  -l-2d1}ei  =  2di-f  dj,  c?  =  (1 +  6)2  +  26di,  d?  =  —d\, 
the  relation  between  the  known  a  and  the  unknown  r  is 


u  -f  v  [ci  +  r(b  4-  r)] 
d  +  t(-c2  +  jt) 


(4.24) 


The  equation  has  a  meaningful  solution  only  if  the  shock  is  attached,  in  which  case  all  the  jump  conditions  at 
the  shock  point  are  evaluated  analytically  and  kept  frozen  during  the  entire  calculation.  The  absolute  value 
of  the  angle  between  the  normal  and  the  upstream  direction  of  the  flow  is  arctan  \/(l  -  r)/r  and  its  sign  is 
opposite  to  the  sign  of  the  deflection.  Since  the  shock  remains  anchored  to  the  corner,  its  velocity  vanishes 
identically.  Therefore,  the  normal  Mach  number  is  Ma,  multiplied  by  che  cosine  of  the  angle  between  the 
impinging  velocity  and  the  normal.  All  the  downstream  values  follow  as  shown  above.  Otherwise,  the  shock 
produced  by  the  wedge  is  detached  and  it  is  evaluated  using  the  same  procedure  as  for  any  other  shock  point, 
with  the  only  constraint  of  the  normal  being  parallel  to  the  wall. 
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III.  SAMPLE  COMPUTATIONAL  RESULTS 


5.  FLOW  OVER  A  WEDGE 

The  code  is  going  to  be  tested  for  cases  of  increasing  complexity,  both  for  inviscid  and  viscous  flows.  In 
principle,  a  single  code  may  be  used  in  all  cases.  Substantial  changes,  however,  must  be  made  to  the  initial 
and  boundary  conditions,  and  it  is  convenient  to  use  separate  codes  when  the  changes  require  more  than  a 
few  IF  statements. 


C  D 


A  E  B 

Fig.  5.1  -  First  wall  geometry 


For  all  cases,  an  orthogonal  grid  is  used.  We  have  used  two  simple  definitions  of  the  wall  geometry: 

First  geometry  (Fig.  5.1)  —  The  lower  boundary  contains  a  rigid  surface  composed  of  three  parts:  a 
straight,  horizontal  segment  between  a  point  E  and  a  point  B,  a  straight  segment  between  point  B  and 
point  C,  and  an  arc  (almost  straight  and  horizontal),  between  point  C  and  point  D.  To  the  left  of  point 
E,  the  boundary  contains  a  straight  segment,  AE,  that  is  considered  permeable  for  viscous  calculations  and 
impermeable  for  inviscid  calculations.  In  the  latter  case,  point  E  has  no  relevance  and  it  is  not  considered 
in  the  code.  The  orthogonal  grid  is  obtained  as  follows.  Let  \p  be  the  wedge  angle,  that  is  the  angle  between 
BC  and  the  horizontal  x-axis.  Inputs  for  the  geometry  are  the  abscissae  of  points  B,  C  and  D  ( z a  being 
assumed  equal  to  zero)  and  the  angle  Let  z  =  x  +  iy  be  the  complex  coordinate  of  the  physical  plane, 
P  =  exp(iV’),  z  =  x  +  iy,  zj  =  Xi  +  iy,, 

Pi  =*/(*- V').  P2  =  7T/(7T+^)  (5.1) 

*1  =[(z  -  xb)/P]p>  (5.2) 

The  function  defined  by  (5.2)  maps  the  broken  line  ABC  onto  a  straight  portion  of  the  x,-axis,  with  the 
origin  at  B.  Let  £  =  (  +  ,'»?  and 

C  =  Wx  -  *ic)P  (5.3) 
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This  function  brings  the  origin  tc  point  C  and  flattens  the  angle  between  BC  and  CD,  leaving  the  segment 
ABC  on  the  (-axis.  We  define 

ox  =  (5-4) 

and  introduce  a  computational  variable,  X,  running  from  0  to  1: 

*  =  «-&)/*.  (5-5) 


Along  the  y-axis,  we  apply  a  stretching,  introducing  a  stretching  parameter,  <ry  and  a  computational  variable 
Y  defined  by: 


■\r  n  ,  tanhffy(Y-l), 

”-yi[1+  tanhtr,  1 

where  Yj  is  a  suitable  value,  to  define  the  thickness  of  the  physical  region.  Consequently, 


(5.6) 


X(  =  1/  ox, 


tanh  ffy 

~  OyYl 


cosh2  [ffy  (Y  -  1)] 


(5.7) 


A  rectangle,  0  <  X  <  1,0  <  Y  <  1,  is  then  mapped  onto  a  domain  that  has  a  lower  boundary  as  shown  in 
Fig.  5.1,  and  the  grid  in  the  physical  plane  is  orthogonal. 

To  have  a  (=constant  grid  line  passing  through  point  B,  after  the  number  of  intervals  along  the  A -axis 
has  been  chosen,  the  value  of  Ax  is  slightly  modified,  without  changing  ax\  consequently,  the  position  of  D 
in  the  physical  plane  is  slightly  modified  from  the  value  initially  chosen.  No  ^constant  grid  line  is  required 
to  pass  through  point  C. 

After  corresponding  values  of  z  and  £  have  been  determined,  or,  G,  <f> i  and  <j>2  are  computed  numerically, 
using  centered  differences. 

The  particular  case  of  a  flat  plate  with  a  leading  edge  at  E  is  trivial. 

Second  geometry  —  A  second  wall  geometry  is  similar  to  the  one  described  above,  but  the  broken  line, 
BCD  is  replaced  by  a  curve,  obtained  by  a  Karman-Trefftz  mapping: 


C-l 

C  +  1 


(5.8) 


after  having  moved  the  y-axis  to  pass  through  D  and  having  rescaled  the  plane  to  let  B  be  at  z  =  -1. 

No  stretching  is  used  on  the  x-direction,  only  a  change  of  scale.  Strong  stretchings  are  used  on  the  y-direction 
to  account  for  the  presence  of  a  boundary  layer.  So  far,  we  have  used  the  stretching  defined  by  (5.6),  with  a 
trial-and-error  preliminary  calculation  to  determine  ffy  after  prescribing  the  height  of  the  first  cells  above  the 
wall  in  the  (-plane  (A??).  A  second  way  of  defining  the  stretching  is  currently  being  tested.  The  object  is  to 
prescribe  a  number  of  intervals  covering  the  boundary  layer,  according  to  an  educated  guess  of  its  maximum 
height  for  a  given  Reynolds  number,  and  to  have  such  intervals  all  evenly  spaced  on  the  (-plane,  using  the 
stretching  only  above  that  strip.  Previous  experience  suggests  a  sensible  reduction  in  computational  time 
being  obtained  by  using  such  a  device. 
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6.  INVISCID  AND  VISCOUS  FLOW  OVER  A  FLAT  PLATE 
Basic  tests  require  calculations  of  inviscid  and  viscous  flows  under  the  simplest  assumptions.  In  any  test, 
residuals  are  defined  as  mean  square  values  of  the  time  derivatives  of  u,  v,  a  and  S,  as  defined  by  (3.4). 
Since  the  code  has  two  levels,  the  second  of  which  corrects  the  derivatives  defined  by  the  first  level,  the  mean 
square  calculation  is  applied  to  the  sum  of  the  time  derivatives  obtained  in  predictor  and  corrector  at  each 
point;  in  this  way,  possible  compensations  of  opposite  contributions  are  accounted  for.  The  logarithms  (base 
10)  of  the  four  residuals  are  plotted  separately.  Note  that  all  values  are  reset  according  to  the  boundary 
conditions  over  the  boundaries  of  the  computational  region,  except  at  the  bottom  line  in  front  of  the  leading 
edge,  where  the  values,  as  updated  in  predictor  and  corrector,  are  accepted  as  valid.  Therefore,  no  special 
care  is  taken  to  obtain  meaningful  derivatives  along  such  boundaries  and  their  values  are  not  accounted  for 
in  evaluating  the  residuals. 

We  begin  with  a  flat  plate.  Obviously,  a  test  of  inviscid  flow  starting  from  an  exact  solution  would 
be  almost  insignificant  because  all  derivatives  would  vanish  identically.  A  test  of  an  inviscid  flow  starting 
from  rest,  instead,  is  quite  significant.  On  the  flat  plate  the  orthogonal  grid  is  obviously  Cartesian.  We 
may,  however,  use  an  unstretched  grid  or  a  highly  stretched  grid.  The  first  test  (Run  343)  was  made  on  an 
unstretched  grid  with  60  intervals  along  £  and  20  intervals  along  rj,  with  the  leading  edge  of  the  plate  at 
£  =  1;  the  length  and  height  of  the  computational  region  are  5  and  2,  respectively.  The  final  =  2.5  is 
reached  after  an  acceleration  phase  lasting  till  to  =  0.5.  “Machine-zero”  values  of  the  residuals  tire  reached 
after  about  500  steps  (Fig.  6.1). 

In  the  second  test  (Run  346),  the  stretching  is  strong  (Arj  =  0.0066).  It  takes  longer  to  bring  the 
residuals  to  machine-zero  (Fig.  6.2)  but  a  flow,  uniform  for  all  practical  purposes,  is  reached  in  about  1000 
steps.  Use  of  a  local  At  (Run  347)  brings  in  a  remarkable  reduction  in  the  number  of  steps  (Fig.  6.3).  Note 
that  during  the  acceleration  phase  a  global  At  is  used,  as  it  appears  from  the  identity  of  residuals  in  Figs. 
6.2  and  6.3  in  the  first  phase  (about  250  steps). 


UIM.IU 


on  ma  xsij  ini  ansa 


Fig.  6  1  -  Residuals  for  run  343 
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Fig.  6.2  -  Residuals  for  run  346 


M1M-3P 


Fig.  0.3  -  Residuals  for  run  347 

For  viscous  flows,  the  most  significant  plots  are  comparisons  of  computed  distribution  of  ti  at  several 
cross-sections  normal  to  the  plate  with  the  theoretical  Blasius  profile.  Since  we  are  dealing  with  compressible 
flows,  we  plot  the  computed  values  of  u/V^  stretching  the  values  of  r\  not  only  by  the  factor  -  Ze) 

(where  £e  is  the  value  of  £  at  the  leading  edge)  but  also  by  the  factor  1/(1  +  64  *  6M^),  according  to 
Stewartson  [2],  The  coefficient,  .64,  was  chosen  according  to  Fig.  3.1  of  [2],  The  comparisons  are  generally 
made  at  three  stations,  called  A,  B  and  C,  respectively,  the  first  is  close  to  the  leading  edge,  the  last  is  close 
to  the  right  computational  boundary  and  the  second  is  somewhere  in  between. 

In  the  following  tests,  the  grid  is  the  same  as  in  inviscid  cases  above.  The  1’randtl  ..amber  is  always 
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equal  to  1  and  the  plate  is  isothermal,  with  the  wall  temperature  equal  to  the  stagnation  temperature  at 
cruising  speed.  The  flow  is  always  started  impulsively.  As  basic  parameters,  we  consider  M «>  and  72.e. 

First,  we  make  m  experiment  to  decide  whether  a  global  At  can  be  used.  For  this,  we  choose  =  0.2, 
He  =  1000  and,  accordingly,  At)  =  0.0066  (run  348).  We  can  see  from  Fig.  6.4  that  the  calculation  converges, 
but  the  rate  of  convergence  is  extremely  low.  Some  differences  in  the  u-distributions  can  be  seen  between 
step  8000  and  step  16000,  not  at  station  A  but  at  stations  B  and  C  (Figs.  6.  through  6.10).  It  seems  that 
the  braking  of  the  flow  produced  by  the  wall  takes  a  long  time  to  propagate  upwards,  due  to  the  fineness 
of  the  grid  in  the  boundary  layer.  Therefore,  the  thicker  the  layer,  the  longer  it  takes  u  to  reach  its  steady 
value.  We  repeat  the  run  using  a  local  At  (Run  349).  The  speed  of  convergence  increases  drastically  (Fig. 
6.11),  and  the  comparison  with  the  Blasius  profile  is  as  good  as  with  the  previous  run  (Figs.  6.12  through 
6.14).  Thus,  we  decide  to  use  local  At  in  all  these  test  runs.  We  also  decide  not  to  attempt  to  fit  shocks,  to 
avoid  introducing  new  stability  problems. 
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Fig.  6.4  -  Residuals  for  run  348 


Many  other  tests  were  made  with  different  Mach  numbers  and  higher  Reynolds  numbers,  varying  Ar; 
and  Y\  occasionally  The  resolution  in  the  boundary  layer  seems  to  be  sufficient,  with  At)  =  0.0066  for 
7 ve  — -  1000,  .001  for  7?.,  =  10000  and  .0005  for  7 le  =  100000.  A  run  made  for  the  highest  Reynolds  number 
and  Ar?  =  001  showed  no  breaking  effects  propagated  from  the  wall  to  the  first  grid  line  above  it! 

The  following  table  shows,  for  ali  the  runs  made  so  far,  an  approximate  value  of  the  order  of  magnitude 
of  the  decrease  in  residuals  in  the  first  6000  steps. 
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Fig.  6.5  -  run  348,  Blasius  comparison  at  A,  step  8000 
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Fig.  6.6  -  run  348,  Blasius  comparison  at  B,  step  8000 


From  such  results,  a  definite  trend  cannot  be  obtained  yet,  although  it  seems  that  supersonic  flows 
converge  faster  than  subsonic  flows  and  runs  with  Ht  =  10000  converge  better  than  others  (perhaps,  because 
the  resolution  in  the  boundary  layer  is  more  appropriate?).  Runs  with  high  Reynolds  number  converge  very 
slowly  and  this  may  hamper  studies  of  more  complicated  cases. 

The  comparison  with  Blasius’  solution  is  generally  very  good.  See,  for  example,  Fig.  6.15  (M„  =  2.5, 
ne  =  10000),  Fig.  6.16  (A/«  =  2,  nt  =  1000),  Fig.  6.17  (A/ro  =  2.5,  Ut  =  100000). 
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Fig.  6.7  -  run  348,  Blasius  comparison  at  C,  step  8000 
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Fig.  6.8  -  run  348,  Blasius  comparison  at  A,  step  16000 
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Fig.  6.9  -  run  348,  Blasius  comparison  at  B,  step  16000 
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Fig.  6.10  -  run  348,  Blasius  comparison  at  C,  step  16000 
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Fig.  6.13  -  run  349,  Blasius  comparison  at  B,  step  8000 
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Fig.  6.14  -  run  349,  Blasius  comparison  at  C,  step  8000 
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Fig.  6.16  -  run  353,  Blasius  comparison  at  C,  step  6000 
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Fig.  6.17  -  run  339,  Blasius  comparison  at  C,  step  4000 
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7.  IN  VISCID  FLOW  OVER  A  WEDGE,  WITH  AN  EXACT  SOLUTION  AS  INITIAL 
CONDITION 

The  first  experiment  in  shock-fitting  computes  an  inviscid  flow  over  a  wedge.  The  initial  conditions  are 
chosen  to  describe  the  exact  solution.  Two  cases  arc  considered,  one  with  an  unstretched  grid,  the  other 
with  a  strongly  stretched  grid  (as  the  one  needl'd  for  viscous  flows).  The  unstretched  grid  is  poorly  suited 
for  a  shock  computation  that  uses  only  “x-type”  shocks  (that  is,  defined  by  intersections  of  the  shocks  with 
^"constant  lines),  if  the  free  stream  Mach  number  and  the  wedge  angle  are  too  high;  that  is  because  the 
shock  tends  to  become  parallel  to  rj  —  constant  lines  and  is  defined  by  too  few  points.  In  this  case,  “y-type” 
shocks  would  be  better  suited. 

To  test  the  code  with  a  stretched  grid,  we  chose  a  case  with  =  6  and  tp  =  25°  (Runs  123  and  124). 
Run  123  (Fig.  7.1)  cuts  the  computational  domain  before  the  upper  corner;  therefore,  no  expansion  occurs. 
Run  124  (Fig.  7.2)  considers  the  entire  domain.  The  results  at  step  2000  are  reasonable  in  both  runs,  but 
in  run  124  there  are  signs  of  uncertainty  at  the  right  boundary,  where  the  shock  reaches  it.  This  is  probably 
due  to  poor  handling  of  the  shock  at  the  right  boundary.  The  rest  of  the  flow  field  seems  to  be  OK  in  both 
cases.  The  residuals  (Figs.  7.3  and  7.4)  reflect  the  inaccuracies  of  run  124;  they  decrease  very  slowly,  with 
oscillations  that  seem  to  increase  in  amplitude.  For  run  123,  the  residuals  reach  machine  zero  in  less  than 
2000  steps. 
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Fig.  7.1  -  Isomachs  at  step  2000  for  run  123 


8.  VISCOUS  FLCW  OVER  A  WEDGE,  WITH  AN  IMPULSIVE  START 
To  test  thr  calculation  of  a  viscous  flow  over  a  wedge,  we  focussed  our  attention  on  the  second  geometry, 
considering  a  flow  defined  by  Mw  =  2.5  and  %t  -  10000.  We  used  a  stretching  defined  by  At?  =  0.001. 
Some  runs  were  made  starting  impulsively  (t0  =  0)  and  running  2000  or  4000  steps  using  local  A t'b  and 
capturing  the  shocks  isentropically  before  introducing  shock  fitting.  Other  runs  were  made  again  starting 
impulsively  and  using  shock  fitting  from  the  beginning.  Interesting  features  were  discovered  during  the 
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7.3  -  Residuals  for  run  123 

debugging  phase.  Shock-capturing  runs  converge  to  machine-zero  residuals  in  less  than  4000  steps  When 
shock-fitting  is  applied,  some  pulsations  appear  in  the  recirculation  zone  and  affect  the  rest  of  the  flow  field, 
interacting  with  the  shocks  that,  in  turn,  oscillate  near  their  origins.  The  upper  parts  of  the  shocks,  instead, 
are  not  affected  by  the  oscillations.  So  far,  the  cause  of  the  oscillations  has  not  been  found  yet;  we  believe 
that,  most  probably,  it  is  purely  numerical  but  a  possible  physical  interpretation  is  not  excluded  Let  us 
make  clear  that  the  calculations  are  still  carried  on  for  not  more  than  10000  steps  and  that  since  we  must 
use  a  global  At  in  the  presence  of  shocks  the  physical  time  span  actually  covered  is  only  of  the  order  of  a 
few  millisecond.  Some  settling  of  minor  oscillations  in  a  different  time  scale  is  not  at  all  unusual  It  is  not 
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Fig.  7.4  -  Residuals  for  run  124 

considered  proper,  however,  to  try  longer  runs  before  making  sure  that  the  code  is  completely  debugged;  the 
logic  of  shock-tracking  and  its  interactions  with  viscous  term  calculations  may  contain  some  hard-to-detect 
traps.  Moreover,  difficulties  may  arise  if  the  computed  region  is  not  sufficiently  high;  shocks  that  reach  the 
upper  boundary  produce  boundary  conditions  that  are  hard  to  implement;  instabilities  are  easily  produced 
at  the  upper  right  corner. 

Here  are  some  results  for  a  typical  calculation.  Fig.  8.1  shows  a  set  of  u=constant  lines  on  the  compu¬ 
tational  plane  at  step  6000.  Note  that  the  V =constant  lines  are  evenly  spaced  on  the  computational  plane. 
Because  of  the  stretching,  the  boundary  layer  occupies  almost  one  third  of  the  region.  Only  lines  defined 
by  ti  _0,  1,  and  2  are  shown  in  the  figure.  The  u=0  line  encloses  a  region  of  negative  u,  which  are  part  of 
the  recirculation  zone.  Two  major  shocks  are  denoted  by  crosses;  the  one  to  the  left  is  rather  weak  and  it 
is  generated  by  the  abrupt  stop  of  the  flow  at  the  leading  edge  of  the  wall;  the  one  at  the  right  is  stronger 
and  it  is  generated  by  the  recompression  at  the  end  of  the  recirculation  region.  The  other  figures  show  some 
of  the  features  in  the  physical  plane.  Fig.  8.2  shows  a  few  u=constant  lines  on  the  entire  computed  region; 
details  of  the  recirculation  zone  are  in  Fig.  8.3;  details  of  the  pressure  distribution  in  the  same  zone  are 

in  Fig.  8.4;  note  the  Constance  of  pressure  on  normals  to  the  wall  and  the  coalescence  of  isobars  into  the 
recompression  shock. 
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Fig.  8.1  -  u=constant  lines  on  computational  plane 
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Fig.  8.2  -  u- constant  lines  and  shocks 
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Fig.  8.3  -  Enlargement  of  Fig.  8.2 
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Fig.  8.4  -  Isobars  in  the  physical  plane 


9.  SHOCK/BOUNDARY  LAYER  INTERACTION 


The  next  flow  considered  was  that  produced  by  shock  reflecting  from  a  boundary  layer.  Only  the  steady 
flow  will  be  considered.  The  free  stream  Mach  number  is  2  and  the  Reynolds  number  is  2.96x  l(P.  A  down  running 
shock  is  produced  by  a  small  3°  deflection.  Figure  9.1  shows  the  isobars  computed  with  the  shock  Fitting  code  and 
fig.  9.2  shows  the  same  computation  with  a  upwind  shock  capturing  code.  The  shock  capturing  code  is  state  of  the 
art  in  that  it  uses  Roc’s  flux  differencing,  for  this  ease  no  limiter  was  required  so  that  the  computation  was  second 
order  everywhere.  A  comparison  of  figs.  9.1  and  9.2  indicates  the  amount  of  shock  spreading  in  the  capturing  result. 
Figures  9.3  and  9.4  show  the  pressure  distributions  at  a  number  of  heights  above  the  plate.  While  on  the  surface  of 
the  plate  the  pressures  look  the  same  off  the  plate  the  shock  spreading  is  obvious.  Figures  9.5  and  9.6  show  details 
in  the  separation  zone.  Figure  9.5  shows  the  long  thin  separation  separation  bubble  and  Fig.  9.6  shows  the  isobars 
and  shocks  in  the  same  region. 
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Fig.  9.1  Isobars  computed  with  shock  fitting  scheme 


Fig.  9.2  Isobars  computed  with  shock  capturing  code 
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Fig.  9.4  Pressure  distribution  computed  with  shock  capturing. 
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Fig.  9.6  Isobars  near  separation 


IV.  CONCLUSIONS  AND  FUTURE  WORK 


The  first  phase  of  our  work  has  demonstrated  the  feasibility  of  performing  accurate  computations  of  shocked 
viscous  flows.  The  computational  scheme  presented  here  is  efficient  enough  for  today's  supercomputers,  and  it  treats 
shocks  in  such  a  way  as  to  avoid  inaccuracies  due  to  pre-  or  post-oscillations,  spreading  and  local  increases  in 
numerical  viscosity.  The  basic  scheme  has  been  demonstrated  on  model  problems  in  order  to  prove  its  soundness. 
Any  future  work  in  this  area  should  continue  to  test  the  scheme  on  flows  which  would  tax  its  capabilities.  In 
particular,  the  flows  considered  thus  far  were  over  somewhat  simple  geometries;  future  work  would  concentrate  on 
more  complex  2-D  flows,  such  as  the  flow  about  a  SCRAM  jet  inlet. 

Only  steady  flows  have  been  considered  thus  far.  The  code  that  has  been  developed  has  a  time  accurate 
capability  built  in  so  that  time-dependent  flows  such  as  inlet  unstart  and  buzz  could  be  considered. 

The  impact  of  fitting  shocks  on  the  accuracies  of  the  computation  of  shock  wave/boundary  layer  interaction 
may  be  the  most  important  aspect  of  any  future  work.  We  have  begun  such  an  investigation  in  the  first  phase  of 
this  work.  An  investigation  of  this  type  may  uncover  problems  with  the  widely  used  shock  capturing  schemes  that 
may  have  broad  implications. 

Only  laminar  flows  have  been  considered  thus  far.  Any  investigation  of  shock/boundary  layer  interaction 
would  be  incomplete  without  the  inclusion  of  turbulence  effects.  These  investigations  would  start  with  simple 
models,  but  they  would  eventually  examine  shock/boundary  layer  interactions  using  complex  turbulence  models. 

The  real  advantage  of  the  scheme  presented  here  may  be  demonstrated  when  three-dimensional  flows  are 
considered.  While  shock  waves  can  be  captured  quite  sharply  (i.e.,  accurately)  in  one  dimension,  they  are  spread 
significantly  in  two-dimensional  flows.  This  problem  becomes  even  worse  in  three  dimensions.  Grid  adaptive 
schemes  can  be  used  some  what  effectively  to  sharpen  captured  shocks  in  2-D,  but  this  has  been  not  been  proven  in 
3-D.  The  complexity  of  three-dimensional  flows  may  require  grid  adaptation  to  obtain  reliable  solutions.  The  grid 
must  be  adapted  to  boundary  layers,  wakes,  vortices,  geometric  irregularities  (comers,  etc.)  and  other  singularities.  If 
the  shocks  are  fit,  then  that  is  one  system  of  singularities  that  does  not  need  grid  adaptation.  This  is  the  appropriate 
system  to  fit  since  shocks  remain  discontinuous  in  the  model  we  are  considering.  The  shock  is  a  few  mean  free 
paths  thick,  which  is  infinitesimal  if  a  continuum  flow  model  is  assumed. 
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